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ABSTRACT 

Numerical hydrodynamics simulations have established that disks which are 
evolved under the condition of local isothermality will fragment into small dense 
clumps due to gravitational instabilities when the Toomre stability parameter Q 
is sufficiently low. Because fragmentation through disk instability has been sug- 
gested as a gas giant planet formation mechanism, it is important to understand 
the physics underlying this process as thoroughly as possible. In this paper, we 
offer analytic arguments for why, at low Q, fragments are most likely to form first 
at the corotation radii of growing spiral modes, and we support these arguments 
with results from 3D hydrodynamics simulations. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities - 
planetary systems: formation — planetary systems: protoplanetary disks 
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Introduction 



The idea that gas giant planets can 



(GIs) in protoplanetary disks (IBos; 



19971 . 



orm rapidly through gravitational instabilities 



19981 ). usually called the "disk instability" theory, 



has now been subjected 



(Pic kett et al 



2002 



2003 



2003, 



1998 



2005 



Rice et al 



; o intensive study with numerical hydrodynamics 



2000 



2007 



2003 



techniques 



2003: 



Gammie 



2005 



Nelson et al. 



2001 



1998 



Ma yer et al 



Mejfa et al. 



2005 



2000 



2002 



Cai et al. 



Nelson 



200- 



2006 



2006: 



2007 



Boss 



2000, 



2001 



Johnson fc Gammie 



Boley et al. 



2006 . 



20071 ). 



Although all investigators agree that massive, cold protoplanetary disks fragment into 
dense clumps under conditions of rapid cooling or local isothermality, there is not universal 
agreement that these fragments are sufficiently long-lived to be considered b ound gas giant 



protoplane t s or t 



2003 



2007 



Boss 



rat fragmentation will occur under realistic disk conditio ns (IDurisen et al. 



2004 



2005 



Rafikov 



2005 



2006 



Pickett &: Durisen 



2007al ). This important 



question may ultimately be decided through simulations alone, but we will gain greater 
confidence in the numerical results if we understand the physics of fragmentation and clump 
longevity through analytic arguments, even if only approximate. 

As a first step in this direction, we consider the special case of a disk evolved under a 
"locally isother mal" assumption, i.e., the disk is assumed to mainta in its initial temperature 



at all positions (IBoss 



1998 



Pickett et al. 



1998 



Nelson et al. 



19981 ) . For brevity, most 



GI researchers refer to this as an "isothermal" disk, even though the temperature of the 
disk is not the same at all positions. Local, thin-disk simulations with a variety of Qs by 
Johnson Sz Gammie ( 2OO3I ) suggest that fragmentation of isothermal disks occurs for Q < 



about 1.4, and the set of all global 3D hydrodynamics simulations for isothermal disks 
referenced in the preceding paragraph generally supports the occurrence of fragmentation 
for such low Qs. The same set of simulations, taken together, also generally confirms that 
disks are subject to growth of nonaxisymmetric GIs for Q < about 1.7 or so and are stable 
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for higher Qs. In simulated disks that exh ibit GIs, the instabilities are in itiated by the 



growth from noise of discrete spiral modes (IPickett et al 



1998 



2000 



isother mal disk simulations, which will be reported in detail elsewhere (IPickett fc Durisen 



2003) . Our own lates t 



2007bl ). show three types of behavior: 1) at low Q, prompt fragmentation of spiral GI 
modes as they first become nonlinear, 2) at intermediate Q, delayed fragmentation in the 
nonlinear regime, probably due to nonlinear interactions of multiple spiral modes, and 3) at 
the highest unstable Qs, nonlinear spirals which do not fragment. 

In this paper, we explore the idea that prompt fragmentation in isothermal disks 
tends to occur near the radius at which a discrete growing spiral mode corotates with the 
disk gas (the "corotation radius" or CR). In Section 2, we give simple analytic arguments 
that compression by spiral shocks must not be too large if the shock-compressed gas is 
going to be able to fragment and that the smallest compression in a spiral mode should 
be found near its CR. In S ection 3, we describe representative results from a simulation in 



Pickett fc Durisen! (j2007bl ) where prompt fragmentation does in fact occur in the vicinity 
of the CR for a discrete fast-growing mode at low Q. The simulation allows us to estimate 
some of the uncertain factors in the analysis of Section 2 and verify that the analytic 



arguments do apply. Section 4 summarizes our main conclusions. 



2. Compression-Induced Stability Against Fragmentation 

Consider cylindrical coordinates (r,<p,z) with the rotation axis of the disk along the 
z-axis, the disk midplane at z = 0, and the sense of rotation in the positive ^-direction. 
In the absence of spiral structures and shocks, a low-Q disk has a half-thickness H <C r in 
the z-direction defined such that Si = 2piH, where Ei is the surface density and p\ is the 
midplane density. We shall assume that the initially fastest-growing discrete spiral mode 
develops more rapidly than any other instability which may give rise to fragments. A spiral 
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shock at r associated with the mode then sweeps material into a sheet of vertical height H 
with thickness L normal to the spiral shock front. Because the gas responds isothermally 



to the shock, there is no vertical jump behind the shock, and the vertical 



remains essentially the same on both sides of the shock (IBoley fc Durisen 



reight of the disk 



2006|). 



We will analyze the gravitational stability of the swept-up gas by assuming that it is 
well approximated by a thin, plane-parallel sheet. This requires L ^ r and 

L < 2a x H, (1) 

where a± is a constant of order unity which accounts for uncertainty about how much 
smaller L must be than 2H for our subsequent analysis to be at least approximately valid. 
Throughout, we will introduce similar factors, each of which is presumably of order unity 
and designated by an atj, where j is an index. They account for uncertainties associated 
with our assumptions. 

If the density p 2 in the swept-up sheet can be considered roughly constant over scales 
of order L, which seems reasonable, then potential energy per unit area of the gas sheet is 

E G = -nGp 2 2 L 3 /Q, (2) 

where G is the gravitational constant. There are several effects that can counter 
the tendency of self-gravity to cause fragmentation, including thermal pressure, tidal 
gravitational stresses, and shearing motions. The Virial Theorem suggests that an infinite, 
isothermal, plane-parallel, ideal gas slab will be stable against fragmentation due to its 
thermal energy alone if \Eq\ < 2Et, where Ex is the total thermal kinetic energy of the 
slab per unit area. This gives 

nGp 2 L 2 < 18a 2 Ci 2 , (3) 
where q is the isothermal sound speed and a 2 accounts for the uncertainty due to deviations 
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from a thin, uniform slab geometry. Stability condition (3) simplifies to 

S s iab 2 /P2 < 18a 2Ci 2 /nG, (4) 

where S slab = p 2 L is the surface mass density of the swept up slab of gas in one of the spiral 
arms. If the spiral mode has m equally-spaced arms that together have swept up a fraction 
/ of the material in each annulus of the disk at the time of fragmentation, then 

Ssiab = 2nrfp 1 /m. (5) 

From equation (4), we see that, for a fixed value of S slab , fragmentation is inhibited 
by large values of p 2 . In other words, strong compression by spiral shocks suppresses 
fragmentation. This somewhat counterintuitive result is peculiar to isothermal slabs. 
Compression by the spiral shock is least near the radius at which the spiral pattern corotates 
with the disk orbital motion. At significant distances from CR, shock compression will 
be strong and will suppress fragmentation. A major aim of this section is to estimate the 
minimum distance away from the CR at which this suppression is effective. 

We first examine criterion (1) by noting that 

L = 2nrfp 1 /mp 2 . (6) 
If the shock is strong and isothermal, then 

P2/P1 = v s h 2 (sin^) 2 /ci 2 , (7) 

where i> s h is the ^-direction speed of the spiral shock front relative to the pre-shock material 
and ip is the angle to the normal of the spiral shock front made by the fluid streamlines 
coming into the shock. At a distance 5r from the CR of the spiral mode, 

^sh = Sa 3 Q5r/2, (8) 
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where 5r/r <C 1 and f2 is the disk's undisturbed angular frequency of rotation. The factor 
as — 1 for a disk in Keplerian rotation but may differ from unity for disks with significant 
self- gravity or pressure support. We define a 4 such that 



H = a^Ci/fl. 

Using (6) through (9), we find that the thin-sheet condition (1) is satisfied if 

2 



(9) 



> 



4nrf 



5r 



(10) 



Now let us consider criterion (3) for co mpression-indu ced suppression of fragmentation. 



For disk gas that is strictly isothermal, the 



Toomrd ( 119641 ) stability parameter Q\ is 



Qi = CiAC/VCrSi, 



where the epicyclic frequency k is \d(r 2 VL) 2 /r^dr] 1 / 2 (jBinney &: Tremaine 



19871) • We set 



(12) 



so that CK5 = 1 for a Keplerian disk. Using (4), (5), (7), (8), (9), (11), and (12) plus the 
definition of Ej, we get that condition (3) for suppression of fragmentation is satisfied if 

2 ' T 2^, f2 



n > 



(13) 



r / 81a2«3 2 «4(sim/>) 2 m 2 Qi 

Together inequalities (11) and (13) show that, for sufficiently large values Sr, both conditions 
(1) and (3) are met. In other words, away from the CR of a growing spiral mode with m 
arms, shock compression due to the mode creates a thin sheet of gas which is stable against 
fragmentation. If prompt fragmentation is to occur, it must do so in the vicinity of the CR. 
According to (11) and (13), as m decreases, the region around the CR where the post-shock 
slab is not stabilized increases. Whether prompt fragmentation will occur thus depends on 
what mode with what corotation radius is likely to reach nonlinear amplitude first. 



- 8- 



So far, there is no rigorous analytic theory which predicts the number of arms for the 
fastest growing nonaxisymmetric mode of a gravitationall y unstable disk. However, a WKB 



analy sis for axisymmetric modes of thin disks with Q < 1 (IToomre 



1964 



Binney &: Tremaine 



19871 ) yields a wavelength for the fastest growing ring-like mode of 



A f w 2.2tt 2 G'Si/k 2 



(14) 



Comparisons between simulations and equation (14) (IDurisen et al.ll2003l ) suggest that, for 
disks unstable to nonaxisymmetric modes, the fastest growing nonaxisymmetric mode has 
a number of arms given by 



rrif = a 6 7ir/\f = a 4 a 5 a 6 Qir/2.2_f/', 



(15) 



where ag ~ 1- Disks do not break up directly into fragments of this size. Instead, the 
simulations show that the mode grows to nonlinear amplitude as a spiral wave. In the case 
of prompt fragmentation, fragments appear first as a smaller-scale instability in the dense 
post-shock region associated with the nonlinear wave. If we set m = rri{ in (10) and (13), 
the thin-sheet criterion (1) becomes 



5r 



> 



8.8nf 



3.1C 



H J 9aia3 2 a4 3 a50!6(snri/>) 2 Qi ctiQi 
and criterion (3) for compression-induced stability against fragmentation becomes 

0.24tt 2 / 2 2.4/C 



(16) 



6r 

I > 



H ) a2a3 2 a4 3 a5ae 2 (smijj) 2 Qf a2a^Qf 



3 ■ 



;i7) 



where 



C = f /a 3 2 a 4 3 a 5 a 6 (smijj)' 2 . 



18) 



The similar numerical quantities on the right hand sides of conditions (16) and (17) 
suggest that, whenever the spiral shock compresses gas into a thin slab, it will be stable 
against fragmentation. However, this also means that the thinness and strong shock 
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assumptions required for use of relations (3) and (7) break down together simultaneously 
near the CR. We cannot invert the inequalities to obtain conditions for instability because 
there is at least one other stabilizing effect, namely, shear. So, strictly speaking, our analysis 
does not tell us whether or when prompt fragmentation does occur, but it implies that, if 
prompt fragmentation occurs, then it must happen in the vicinity of the CR for conditions 
under which the inequalities are reversed. We can see from the Qi-dependence of (17) that 
this is more likely to happen for low Q ; . 



3. Comparison with Simulations 



3.1. Methods and Initial Conditions 



Adopting 



Pickett et al. 



the s ame basic star/disk model and 3D hydrodynamics code used by 



(119981 ). we have computed a series of locally isothermal simulati ons of disks 



with d ifferent, but constant, values of Q. The r, z-resolution is the same as in 



Pickett et al 



(|l998l ). but the ^-direction resolution is increased from 128 to 512 azimuthal cells in 27r 



radians. T 
satisfy the 



lis seems sufficient to reso 



Truelove et al. 



fll997h and 



ve frag mentation in modes of moderate m and to 



Nelson! (120061 ) criteria for avoiding purely numerical 



fragmentation prior to the onset of fragmentation. So far, the series of simulations includes 
Q = 1, 1.15, 1.25, 1.35, 1.5, 1.6, and 1.7. For reasons explained in Section 3.2 below, Q 
here, without the subscript "i", is computed using (11) but replacing the isothermal sound 
speed by the adiabatic sound speed for an ideal gas with ratio of specific heats 7 = 5/3. 
Consequently, 

Q = (5/3) 1/2 Qi = 1.29Qi (19) 
The star /disk model is the same massive "stubby" disk used in 



Pickett et al. 



(1998. 



2000), with disk-to-star mass and radius ratios of about Md/M s = 0.25 and Rd/R s 
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7.1, respectively. In this case, however, in order to avoid the severe cost of the Courant 
time-step limitation near the rotatio n axis, the disk is de tached from the star using the 
localized cooling method described in 



Pickett et al. 



(120031 ). and the stellar mass distribution 



and gravitational potential are frozen. 



3.2. Prompt Fragmentation at Corotation 



In the full Q-survey, which will be reported in greater detail elsewhere 



(IPickett fc Durisen 



2007bl ). the initial disk is given a small amplitude, random, cell- 
to-cell density perturbation. This allows fast growing modes to organize themselves and 
grow in a few dynamic times. The three behaviors described in the Introduction are 
found over the following ranges: 1) 1.0 < Q < 1.25: Prompt Fragmentation. Fragments 
appear as soon as growing modes reach nonlinear amplitude. 2) 1.35 < Q < 1.6: Delayed 
Fragmentation. Discrete modes do not fragment as soon as they reach nonlinear amplitude; 
instead, the appearance of dense fragments is delayed by several to many pattern rotations 
and appears to be associated with nonlinear interactions of different patterns or arms. 3) 
Q = 1.7: Nonfragmenting Spiral Arms. The disk is unstable and develops strong spiral arm 
structure, but no fragmentation occurs over the duration of a long simulation. We have 
not yet tested the upper bound of GI stability for this star/disk model, but we do know 



fromlPickett et al. 



(119981 ) that a Q = 2.0 disk is stable against growth of nonaxisymmetric 



structure, so the Q-limit for onset of GIs is somewhere between 1.7 and 2.0. Clearly, the Q 



below which fragmentation occurs for these isothermally evol ved disks is between 



Johnson fc Gammid (120031 ) in 



.6 and 



1.7. We compare this fragmentation limit with that found by 
Section 3.2 below. 

For the low Qs of interest here, where prompt fragmentation occurs, many modes grow 
rapidly at once, and it is difficult to discern which mode or pattern period is associated 
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with a give n fragment or set o f 



techniques (IPickett et al 



1998 



fragm ents. So, for this paper, we apply Fourier analysis 



20001 ) to the Q = 1.25 case and determine, as best as we 
can, that the fastest-growing mode is a particular m = 4 mode. We then run an additional 
simulation with a pure, low-amplitude Sp/p ~ cos40 density perturbation over the radial 
range where the fastest growing m = 4 pattern was detected in the random perturbation 
simulation. As a cautionary note, we point out that it is difficult to be sure we have truly 
isolated the most unstable mode, because many modes grow with similarly fast growth 
rates. A midplane greyscale of the disk close to the moment of prompt fragmentation in the 
run with the pure cos40 hit is shown in Fig. 1, where it is apparent that the fragments do 
indeed appear near the CR of the growing mode. 

We can use the simulation to estimate some, but not all of the various parameters that 
enter into our analysis in Section 2. The rotational shear of this massive disk is somewhat 
non-Keplerian, so that 0:3 ~ 0.9 and a§ ~ 0.9. Using Si and the midplane values of pi, C;, 
and Q near the CR in (9), we find that ~ 0.6. Use of these Oj-s and the known values 
of Qi, H, the r of CR, and mj = 4 in (15) gives a§ ~ 0.7. It is difficult to measure ip in 
our Eulerian code, but it is relatively easy to determine the pitch angle of the spirals to be 
about 20 degrees prior to the onset of fragmentation. So we use ip m 20 degrees as a crude 
estimate, although the pre-shock ip is likely to be somewhat larger than the pitch angle due 
to radial motion of the fluid as the wave becomes nonlinear. The fraction of mass / in the 
spirals just prior to fragmentation is also difficult to determine precisely but appears to be 
about 0.2. Using these values in (18) gives ( ~ 16. If we assume «2 = 1, then (17) becomes 

5rV 11 



HJ>Q!' < 20 > 

Unfortunately, it is difficult to determine the true 5r for the parts of the spirals that go 
into the dense fragments, and so it is difficult to verify that (20) is precisely satisfied. 
Nevertheless, for Q = 1.25 (Q\ = 0.97), (20) gives 5r ~ 3 to AH for the edge of the stabilized 
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region, or about six cell widths from the CR. The clumps are well within this distance from 
the CR in disk midplane images near the time of fragmentation. 

Other aspects of our analysis can also be checked against the simulation results. For 
instance, the thinness of the arms in the midplane view of Fig. 1 suggests that condition (1) 
is satisfied, but this observation does not give us much insight into the proper value of a.\. 
Finally, consider mj. Both Qi and H in (15) vary linearly with q, while other parameters 
involved in Qi and H, like k, Q, and Si, are determined primarily by the disk's radial 
structure, which does not vary much with c\ for cool disks. So we expect rrif to be relatively 
insensitive to Q until the disk becomes hot enough for radial pressure gradients to play 
some dynamic role. In fact, our set of simulations with v aried Q indicate that m = 4 is 



2007bl). The 



the fastest growing mode for all Q between 1.15 and 1.5 (iPickett &: Durisenl 
Q = 1.0 case has so many fast growing modes with various ms, including m = 4, that we 
cannot determine which of them grows most rapidly. 



3.3. Fragmentation Criteria 



Comparison of fragmentation criteria derived by different authors involves some 



discussion about how Q is evaluatec 



adiabatic sound speed (IPickett et al. 



For our disk models, we use (11) but with the 



1998 



2000 



20031 ) . One of the goals of our body 



of work has been to consider the same basic equilibrium disk models evolved under 
different assumptions about the equation of state of perturbed fluid elements. Because 
the models are derived from isentropic equilibrium configurations, it makes sense to use 
the adiabatic sound speed to compute a common reference Q that uniquely designates 
individual disk models within a set of related models. For the purposes of comparison with 



other work, h owever, we need to c o mput e Q t 



that we can. 



Johnson fc Gammid (120031 ) and 



le same way other authors do, to the extent 



Mayer et al. 



(12004] ) use what we call here 



13 



Mayer et al 



Q\. I Johnson &: Gammid (120031 ) report fragmentation if and only if Qi < about 1.4; the 



(120041 ) paper suggests a somewhat larger value. We find fragmentation setting 
in for Q < 1.7 or, equivalently, for Q[ < about 1.32, which is in reasonable agreement 
with these other works, especially considering that the precise limit in global simulations 
probably depends somewhat on the detailed structure of the disk. 

Dynamically, it may seem obvious that the isothermal sound speed is the appropriate 
one to use in this case. However, the perturbations in our simulations are locally isothermal, 
i.e., the temperatures are fixed spatially throughout the grid. Sound waves traveling in 
the (f) direction are truly isothermal, but waves moving in the r or z directions are not. 
Thus, it is not absolutely clear what sound speed should be used in our approximate 
stability relations, nor how truly "isothermal" the spiral shocks will be. This uncertainty 
is not explicitly accounted for by any of the <x,-s in our analysis. It also adds an unknown 
level of uncertainty when comparing gird-based locally isothermal simulations with SPH 
simulations in which the temperatures of the particles are kept fixed instead and thus may 
more accurately represent fluid elements that are evolving isothermally. 



4. Conclusions 

We have used simple anayltic arguments to estimate where fragments can first appear 
in a low-Q gravitationally unstable disk when the spiral shocks are locally isothermal. 
Assuming that the disk behavior is dominated by a single, coherent spiral mode and that 
the dense gas behind the spiral shock can be well approximated as a thin slab behind 
the shock, we have shown that compression from isothermal shocks will stabilize the slab 
against fragmentation away from the CR. The minimum distance from the CR at which 
the stabilization becomes effective gets larger for smaller Q[. So the analysis suggests 
that, if fragments are going to form promptly upon nonlinear growth of a mode, then 
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this will happen first near the CR of the mode and only for low enough Q\. The result is 
prompt formation of high density blobs near corotation in each spiral arm of the dominant 
pattern as it becomes nonlinear. A locally isothermal 3D hydrodynamic simulation of 
a well-studied star/disk model with Q = 1.25 (Q\ = 0.97) confirms this prediction. We 
plan to attempt similar analyses in the future to understand fragmentation under different 
conditions, including delayed fragmentation and fragmentation in radiatively cooled disks. 
We suspect that, as in the isothermal case, prompt fragmentation in radiatively cooled 
disks will tend to occur under more extreme conditions tha n fragmentation itself. In other 



words, if fragmentation generally occurs when t coo iQ < c 7 / (IGammie 



2005 



Mejia et al. 



2001 



Rice et al. 



2003 



20051 ) . where c 7 / is a constant that depends on the adiabatic index 7 and 
where t coo i is the time it takes the disk to radiate away its internal energy, then we would 
expect prompt fragmentation to occur when t coo ifl < c 7P , where c 7P < c 7 /. 

point we have made elsewhere 



Pickett et al. 


, — 
1998. 


200(J. 


2003; 



Durisen et al 



2003 



Mejfa et al. 



2005 



Pickett fc Durisen 



2007al ). The occurrence of fragmentation in disks does not in itself demonstrate that 



protoplanet formation by disk instability is really possible. The key issues are whether 



conditions for dis 



2006 



t fragmenta tion are achieved in real disks (IRafikov 



Boley et al 



2005 



2006 



Cai et al. 



2006 



20071 ) and whether de nse clumps, if they do form , are able to evolve 



into protoplanets over many orbital periods (IPickett fc Durisen 



2007al ). Although our 



criteria suggest where potential precursors to gas giant planets may appear under certain 
restrictive conditions, they do not say anything about the longevity of the clumps once 
formed. 
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Figure Captions 



Figure 1. Midplane mass density grayscale. Shown is a grayscale intrepretation of the 
midplane mass density at the time when fragments first appear in the 3D hydrodynamics 
simulation of a Q — 1.25 disk given a cos40 initial perturbation. The greyscale spans six 
orders of magnitude in density. The circle indicates the location of the corotation radius 
(CR). Note that the first clumps to appear are quite close to the CR for the four-armed 
mode. 



